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Abstract 

We present a two-dimensional granular model for the mechanical behavior of an 
ensemble of globular grains, during solidification. The grain structure is produced by 
a Voronoi tessellation based on an array of predefined nuclei. We consider the fiuid 
fiow caused by grain movement and solidification shrinkage in the network of channels 
that is formed by the faces of the grains in the tessellation. We develop the governing 
equations for the fiow rate and pressure drop across each channel when the grains are 
allowed to move, and we then assemble the equations into a global expression that 
conserves mass and force in the system. We show that the formulation is consistent 
with dissipative formulations of non-equilibrium thermodynamics. Several example 
problems are presented to illustrate the effect of tensile strains and the availability of 
liquid to feed the deforming microstructure. For solid fractions below Qs = 0.97, we 
find that the fiuid is able to feed the deformation at low strain, even if external feeding 
is not permitted. For solid fractions above Qs = 0.97, clusters of grains with "dry" 
boundaries form, and fiuid fiow becomes highly localized. 
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1 Introduction 



The last stage solidification of alloys is a critical step in casting and welding processes during 
which several defects can form [1] . In addition to porosity, the most severe of these defects is 
probably hot cracking, i.e., a spontaneous failure of the material while it is still semi-solid. 
This defect, which typically occurs in dilute and hard alloys, limits the productivity of cast 
houses and restricts the range of alloys that can be produced. [2] It also severely limits the 
weldabihty of this class of alloys. 

In dilute alloys, solid grains are separated by thin continuous liquid films up to high vol- 
ume fraction of solid gg (typically up to gg ~ 0.95), especially at high angle grain boundaries, 
where coalescence of solid grains is made difficult by the large grain boundary energy [3] . As 
stresses and strains are generated in the not yet fully coherent sohd, the presence of these 
liquid films together with the low permeability of the mushy zone near gs — 0.95, which 
prevents efficient feeding, makes the material extremely brittle. [4, 5] Deformation induced 
by thermal and solidification shrinkage tends to localize at these liquid films, which act as 
a brittle phase. They pose little resistance to tensile strains and, as they open, the newly 
created volume cannot be fed by intergranular liquid. The so called "brittle range," meaning 
the Qs-iange where the mushy zone exhibits low strength and low ductility, can be measured, 
for example, by tensile tests on partially solidified alloys. [6-8] 

Accurate prediction of hot cracking requires the knowledge of: 

1. The thermal history of the cast or welded part. This is usually done using volume- 
averaged multi-phase heat and mass transfer models [9, 10]; 

2. Realistic mechanical constitutive equations for the semi-solid alloy [7]; and 

3. Conditions under which a hot tear will form, also known as a hot cracking criterion 
[11-13]. 

Hot tearing models based on such data have been implemented successfully in casting sim- 
ulation codes by several authors. [14-16] 

These existing models treat the material as a continuum whose properties are represented 
by volume-averaged quantities at a length scale that is large compared to the microstructure, 
and small compared to the variations of macroscopic fields. By their very construction, they 
cannot include the important transitions in microstructure at high volume fraction solid, 
such as the localization of liquid films to the periphery of increasingly large grain clusters, 
as described by percolation theory. [17] Thus, the loss of locality of the length scale implies 
that no appropriate representative volume element can be defined. This limits the utility 
of such models, because it is in precisely this transition region that hot cracks form, as 
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recently demonstrated experimentally by Gourlay and Dahle, emphasizing the importance 
of the granular nature of solidifying alloys. [18, 19] 

In order to go beyond the limits imposed by the volume-averaged methods, in this work 
we adapt granular or discrete element models (DEM), which have been developed in the 
context of the mechanics of granular materials. [20-22] These models simulate the behavior 
of a large number of spherical grains, which may be either rigid [20,21] or deformable [22]. 
They consider the interactions between the grains due to solid-sohd contact but neglect the 
influence of the surrounding medium. To extend these models to the case of semi-solid alloys 
at high solid fraction, we must consider not only the solidification of each grain, but also the 
infiuence of fiuid fiow on the mechanical behavior of the mushy zone. 

The first model using the granular approach for the solidification and coalescence of 
globular grains in 2D was proposed by Mathier et al. [23] This work was further developed 
by two of the present authors. [24] In [25], a percolation analysis was presented that identified 
the various transitions in the mushy zone that appear naturally in this approach, as well as a 
model for liquid feeding. However, all these contributions made the major assumption that 
the grains remained fixed. In the present contribution, we derive a 2D granular mechanical 
model for globular microstructure that includes both the flow of intergranular liquid and the 
displacement of solid grains. 

2 Solidification model 

The 2D mechanical model that we present here is based on the solidification model presented 
in detail in previous publications. [24,26] The model is appropriate for inoculated alloys whose 
final grain structure is fine and globular. For the sake of completeness, we briefly recall the 
main features of this sohdiflcation model before presenting the mechanical DEM approach. 

Consider a population of grains distributed randomly in space with a specifled average 
density. Nucleation is assumed to be instantaneous, and the temperature difference across a 
typical grain is taken to be small compared to the undercooling, i.e., small thermal gradient. 
The flnal grain structure is therefore close to the Voronoi tessellation of the original set of 
nuclei. [27] Beginning with a random distribution of nucleation centers, the Voronoi tessel- 
lation of the set is computed from the mediatrix of the segments connecting two neighbor 
nucleation centers (see Fig. 1). A 2D grain is then identified as an ensemble of triangles 
having the nucleation center as one corner and two vertices on the grain boundary. Com- 
putation in the model proceeds in two steps. In the first step, the interface of each grain 
during growth is approximated by a linear segment parallel to the grain boundary. In this 
step, we neglect solute fiux between elementary triangles, which reduces the solute diffusion 
calculation to a one-dimensional problem in each triangle. We assume complete mixing in 
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Figure 1: Left, grain as computed with the Voronoi tessellation method. Right, notations 
for a liquid channel used as the basic element of the mechanical model. 

the liquid. Back-diffusion in the solid can be easily incorporated in the model. [24,26]. In 
the second step, the sharp corners of the polyhedral grains are smoothed using a procedure 
that accounts for the local curvature undercooling. Further details of this model, together 
with a discussion of its limitations and domain of validity are given in Ref . [24] . 

3 Mechanical model 
3.1 Basic hypotheses 

We begin the mechanical model with the microstructure computed by the solidification 
model. The idea is to derive a DEM for the mechanical response of the mushy zone that 
retains the main physical aspects at the scale of the grains, yet is simple enough to be 
computationally tractable. To that end, the following assumptions are made: 

1. The grains are undeformable. This assumption is adequate to obtain insight into the 
physics of hot tearing, where deformations are small. A more complete model of the 
rheology of the mushy zone should also include deformation of the solid, which is known 
to be important for volume fractions of solid Qs ^ 0.6, known as the traction coherency 
sohd fraction. [7,28,29]. 

2. Grain movement can occur only by translation. This hypothesis seems very restrictive, 
but experimental studies have shown that at high temperature and high solid fraction, 
the most important deformation mechanism is grain boundary sliding. [6] Making this 
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assumption greatly improves the computational efficiency of the model, because detec- 
tion of contacts between polygonal grains is difficult if rotations are permitted. [20] 

3. Liquid channels smaller than a pre-defined coalescence interaction distance S, on the 
order of a few nm, are taken to he fully solid. [3,23]. Using such a cut-off improves the 
computational behavior of the model by eliminating very large coefficients that would 
otherwise appear in the matrix of the linear system derived in Sect. 3.5. 

4. The intergranular fluid is Newtonian and incompressible, and no-slip conditions apply 
at the solid-liquid (s — £) interface. We adopt this constitutive model for the fluid 
even though the channels between grain arc very thin. We note that measurements by 
Israelachvili showed that the viscosity observed in bulk samples is still valid for films 
as thin as 5 nm. [30,31]. Another effect observed in thin films is slipping at s — £ 
interfaces which may be due to the formation of a nanometer scale air gap. [31,32] 
This phenomenon occurs when the shear stress exceeds a critical value, and might be 
an interesting phenomenon to consider in hot tearing, especially if gas porosity appears 
at grain boundaries. It is not included in this work. 

3.2 Notation 

Our model considers the liquid network formed between the grains. The basic element is a 
liquid channel surrounded by two solid grains, designated a and b in Fig. 1. The element has 
four integration points. The first two points arc the grain centers, denoted Oa and Ob for 
grains a and b, respectively. At these points, we consider two conjugate vector quantities, 
the velocity of the grain {V}q^ (or {V}q^) and the force exerted by the grain on a liquid 
channel {F}^ (or {F}^). The two other integration points are the ends of the liquid channel, 
denoted i and j. The conjugate quantities considered at these points are the fluid flux 
(or and the pressure Pj (or Pj). Uppercase letters are used for these entities at the 
integration points, whereas lowercase letters will be used for the associated fields (e.g, {v} 
for the velocity and p for the pressure at any point in the liquid channel) . 

Relative translation between neighboring grains can produce a mismatch at the extremi- 
ties of the channel. With reference to Fig. 1, we introduce the following length measures: L 
is the length of the s — £ interface for one channel, Lc is the length of the channel where the 
two grains effectively face each other, and Lja is the length of the s — £ interface from the 
center of the channel to the extremity of grain a near vertex j. We also define Lj^ — L — Lja- 
The half- width of the channel is designated as h. Note also that even though the channel 
can be curved at its extremities, it is modeled with straight lines for the calculation of the 
pressure field. ^ The special case where relative motion of the grains causes their faces to no 
^Note that the fltix balances that we introduce in the next section with the polyhedral envelope of the 
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longer align anywhere, i.e., Lc — 0, is treated in Sec. 3.6. 

The volume change associated with solidification shrinkage, i.e., due to the density dif- 
ference {ps — pi) between the solid and liquid phases, produces a compensating fiow in the 
liquid near the interface. It is most convenient to develop the expressions for the flow in 
terms of the normal and tangential components of the velocity, denoted with superscripts 
"n" and "f , respectively 

{VL=f"'0 = f 'o- \ (1) 

° \V: ) \VS.-I3V) 

where (3 is the sohdiflcation shrinkage (/3 = PsI Pi — 1), and v* is the speed of the s — I 
interface as given by the sohdiflcation model. We choose the element normal vector {n} to 
point toward the liquid from grain a, and the tangential vector {t} is directed from vertex i 
to vertex j. Note that, due to this orientation convention, the fluid velocity at the interface 
with grain h is written: 

Symmetry of the sohdiflcation model implies that the velocity v* of the s — ^ interface of grain 
h is equal to that of its neighbor grain a. The coordinates in the {t} and {n} directions are 
noted X and Y, respectively. The origin is deflned at the center of the channel (see Fig. 1). 

3.3 Integration of the constitutive equations 

The scaling analysis given in Appendix A shows that the X— direction momentum balance 
in the channel reduces to the following simpler relation between the pressure in the liquid 
channel p{X) and the fluid velocity 

dp _ d\x 

where p, is the dynamic viscosity. The volumetric flow rate in the channel (j)i^j{X) from 
vertex i to j is deflned by 

<Pi^j{X)= [ vx{X,Y)dY (4) 

J-h 

We also have the continuity equation for a constant density fluid in 2-D, 

dvx , dvY . . 

dX+W^^ 

Integrating Eq. 5 over the width of the channel and using the deflnition of in Eq. (4) 
gives 



grains are strictly the same as those done for the rounded grains. Indeed, as the fluid will be considered as 
incompressible, the mass of fluid in between the grains and their polyhedral envelope remains constant. 
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where 



a ^ Oh 



(7) 



represents a source or sink in the channeL Note that if the velocity of the grains is zero 
(or uniform), i.e., V'' = 2Pv*, Eqs. (3) and (6) reduce to those given in [24], where we 
considered only the shrinkage-induced intergranular flow. 

Equations (3) and (6) are integrated in Appendix B for the general case. To provide 
a better understanding of the underlying physics, let us consider a few special cases where 
the grain and/or liquid movement isolates the individual contributions of certain important 
phenomena. Since the governing equations are linear, solutions for cases that combine these 
phenomena can be obtained by superposition of these simple cases. The cases are illustrated 
in Fig. 2. 



Constraints at 
integration points 



(b) 



(c) 



V" 



Resulting actions 



Due to 
constraints difference 



F?=h(f]-Pi) 



$.=2h_(p-p) , 



^1 ibb laa/ 



/^J jy'b "-ja'a 



I^Fb" 







Due to 
constraints average 



5 =-^^ ^^ 



F" =L. P+L. R 

a la I ja J 




Figure 2: Representation of the various individual constraints imposed on a liquid channel 
and on the grains: (a) Simple flow in between two immobile grains; (b) Normal displacement 
of the grains; (c) Tangential displacement of the grains. 



Case (a), corresponding to the top hne of Fig. 2, neglects shrinkage- induced flow {ps — pi), 
and considers an imposed pressure differential between vertices i and j. The grain velocities 
are set to zero for this case. In the part of the channel where the two s — i interfaces face 
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each other, the fluid velocity profile vx{Y) has the usual parabolic profile given by 

Thus, the volumetric flow rate at the vertices i and j is given by: 

^^^^^'^^■^ = 3iJl,^^' ~ ^ -^.■(^.,^.) (9) 

This equation describes the flow between two parallel planes. In 2D, any contact between 
neighboring grains blocks all flow in the channel, whereas in 3D the liquid can flow around 
the contact. To include this important 3D feature in our 2D model, we assume that when two 
grains make contact, a small "pipe" remains open. The radius of this pipe, r^, is estimated 
to be of the same order of magnitude as the radius of curvature at the grain corner derived 
in [24] (see Fig. 1). Moreover, if the solidiflcation microstructure were to be extended in the 
third dimension, there would be one such pipe each average grain diameter Dav Summing 
the flow rates for the channel and the pipe gives 



(10) 



where we have represented the flow rate in the pipe using the standard Poiseuille solution. 

We show in Appendix A that the pressure loss in the part of the channel where the 
s — £ interfaces of the two grains do not face each other can be neglected in comparison to 
the remaining terms. Therefore, the pressure is constant in these locations, and decreases 
linearly from to Pj along Lc- This simple pressure profile in the liquid channel can be 
integrated to obtain the forces exerted by the grains on the liquid to yield 

{F}.(P..P,)=( ) (11) 

and 

{FK(P.,P,)=( l^^-f ) (12) 

Of course, an opposite force is exerted by the fluid on the grains. 

Next, we consider the individual effect of a grain displacement in the normal direction, 
with Vq^ — Vq^ — and Pj = P,- = (Fig. 2, case (b)). Equations (3) and (4), combined 
with a no-slip condition at the s — £ interfaces, link the pressure gradient in the X— direction 
to the fluid flow. 



With the help of Eq. (6), we have then: 

2h^ 



3/x dX 



p V,- - v: (14) 



Integrating Eq. (14) twice in X, and considering the imposed symmetry on the pressure gives 

3/^ - Va) 



2h 



(15) 



This term represents the change in pressure induced by the fluid flow required to compensate 
the channel expansion {V^^ > V^) or constriction {V^^ < V^). Equation (15) can be integrated 
once more over the length of the channel to obtain the forces exerted by the grains on the 
liquid: 



-^i (Lj2hf (n" - K") 



{F}aiv:,vn = { ,.,,3.... .... I (16) 



and 

{^>''''«"'^"'^(.(W2/.m"-K») I '''' 

The symmetry of the formulation implies that at the center of the channel 

0,_,i^^o(Kr,n") = o (18) 

As the fluid is incompressible, a flux balance is readily written at the channel vertices as 

UV:,V,^) = L,,V,"-L,aV: (19) 

<^>,{v:,vn = L,,v,'^-L,av: (20) 

The final case that we consider isolates the effect of tangential displacement of the grains 
(Fig. 2, case(c)). This situation corresponds to pure shear of the liquid channel, and thus 
we have: 

{FLK.n')=(-''<^«-^"M (21) 



and 



Note that we have neglected the shear forces in the portions of the channel where the two 
grains do not face each other. Considering the flow relative to the initial (reference) conflgu- 
ration, the average tangential displacement of the grains induces a fluid flow at each vertex 
given by 

MV:,V,')^2h^^^±^ (23) 

^MX)--2hy^^ (24) 

If — —Vij, there is no net flow in the channel, i.e., the liquid experiences a perfectly 
symmetric shear stress. 
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For a more general situation in which all three of the phenomena just discussed may 
occur, the various contributions to the fluid flow and to the forces on the grains can be 
written as the sum of the individual simple cases just described, with the result 

IF} = I ~ ^'^^ /^(^c/2/i)(n* - K*) ^ .26) 
^i={^ + ^7^] (Pi - Pj) + ^h^t^ + L,,n" - L,„C (27) 



4 \ T/t I T/t 



We now proceed to develop these expressions into a form suitable for a finite element formu- 
lation. 



3.4 Elementary matrix 

In order to obtain a matrix form, we collect the primitive variables into a vector {U}"^ = 
{Pi, Pj, V^, Vj", V^, Vff ), and the associated flow rates and forces into a second vector {W}^ = 
$j, F^, F^, F^, Fl). The superscript "T" indicates the transpose of the vector. Equations 
(25)- (28) can then be written in matrix form as 



{W} 



where 



a 



[E] {U} 





+Ci 


-Ci 


~Lia 


+Lib 


+h 


+h 


\ 


/ 


Pi 


\ 




-Ci 


+Ci 


~Lja 


+Ljb 


-h 


-h 






p^ 








-\-Lja 


+C2 


-C2 


















—Lib 




-C2 


+C2 


















-h 


+h 








+C3 


-C3 










v 


-h 


+h 








-C3 


+C3 


/ 


\ 




/ 



2/i3 



+ 



3/iLc S^DavLc 



C2 = /X 



2h 



2h 



(29) 
(30) 



Please note that, unlike the usual formulation of mechanical problems using the finite clement 
method, we have chosen to group together the normal and tangential components of the forces 
and velocities acting on grains a and b. One should also keep in mind that the component 
= Vq^ — /3v* (or V^^ — Vq^ + Pv*) includes both grain displacement and sohdification 
shrinkage. 

It is interesting to compute the power dissipation fl in the channel 



^ = {V}, ■ {F}, + {V}„ ■ {F}„ + P,$^ 



(31) 
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or equivalently by: 

Q = {U}^ [E] {U} 

Expanding Eq. (33) using Eq. (29), we obtain 
2h' 



2h 



n\2 



t\2 



(32) 



(33) 



Thus, the element matrix of the system is positive definite and represents a quadratic form 
related to dissipation in the channel. 

It is convenient to decompose the matrix [E] as the sum of its symmetric and antisym- 
metric parts, [S] and [A], respectively 



[S] 



and 



+Ci 


-Ci 











-Ci 


+Ci 

















+C2 —C2 














—C2 +C2 














+C3 


-C3 








0- 


-C3 


+C3 










+h 


+h 








~Lja '\'Ljb 


-h 


-h 




+Lja 











—Lib 













-h 


+h 











-h 


+h 












(34) 



[A] 



V 

Inserting this decomposition into the power dissipation yields 

Q = {Ur([S] + [A]){U} = {Ur[S] {U} 



(35) 



(36) 



since {U}^ [A] {U} = 0. The antisymmetric part of [E] , [A] , represents coupling between 
the pressure in the fiuid channels and the displacement of the grains. This term does not 
dissipate energy in our formulation. 

This form of the matrix is consistent with the Onsager-Casimir theory of transport phe- 
nomena. [33] Consider the two sets of conjugate quantities. One set is invariant to time 
reversal (Pj, P,-, {F}^ , {F}^), whereas the other set (^i, ^j, {V}^ , {V}^) changes sign with 
time reversal. The symmetric matrix [S] couples quantities that have different behavior with 
respect to time reversal. The anti-symmetric matrix [A] couples the quantities with the same 
behavior. For example, is related to Pi,Pj by [S] and is related to {V}^,{V}^ by [A]. It is 
interesting to note that this fundamental relation is obtained by a simple integration of fluid 
flow equation [34]. The usual volume- averaged formulation with the same set of unknowns 
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for describing semi-solid materials (pressure in the liquid, velocity or deformation rate in the 
solid) does not reveal the symmetries inherent in the present model. 

3.5 Global problem 




Figure 3: Construction of the global problem from the elementary matrices: The sum of 
forces on each grain is zero and the sum of fluxes flowing in and out of each vertex is zero. 



To assemble the global problem from the contributions of the individual elements, both 
the unknown velocities at the grain centers and the resulting forces are expressed in a global 
coordinate system {X,Y) (see Fig. 3). One has for example: 



(37) 



where Nx and Ny are the components of the channel normal vector in the global frame. 
The pressures and fluxes for each channel in the global coordinate system are obtained by 
the usual transformation rules, i.e., 





{Wy = [Q] {W} = [Q] [E] {U} = [Q] [E] [Q]^ {U}' = [E]' {U}' 



(38) 



where ({WjT = , F,^ , , F^), {{UYf = , n"", KT, H^) and [E]' = 

[Q] [E] [Q]"^. The transformation matrix [Q] is given by: 



[Q]- 



/ 


1 
























1 
























nx 





—ny 


















nx 





—ny 












ny 





nx 







V 











ny 





nx 


J 



(39) 
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Since T/" = V^^ - i5v* and = Vg^ + Pv* , Eq. (29) in the {X,Y) frame becomes 



<T). 

P^ 



= [Ql |E] IQl 



T 



( P. \ 



P3 



[Q] [E] 



V 






+pv* 





(40) 



The last term, which wc caU {B}', is associated with sohdification shrinkage and is known, 
coming from the external solidification model. The balance of mass and force is obtained 
by summing Eq. (40) over all elements and setting the result to zero. This procedure is 
similar to the matrix assembly in the standard finite element method, i.e., each contribution 
associated with a given liquid channel with a local numbering (a, b, i, j) (Sec. 3.4) is added 
to the global matricial problem with a global numbering of all the grains and vertices. The 
result is written in the compact form 



[Etot]'{Utot}' = -{Btot}' 



(41) 



where the vector of unknowns {Utot} contains the velocities of the Ng grains (V^^, VqJ in 
the global frame and the pressures Pi at the Ny vertices. 

Integration points located at the external boundary of the global domain are subject to 
boundary conditions. An imposed flux on a channel vertex and an imposed force on a grain 
are Neumann-type boundary conditions in this formulation. These are taken into account 
by adding the imposed constraint into the global vector — {Btot}'- Boundary conditions 
specifying either the velocity of a grain or an imposed pressure at a channel vertex are 
essential boundary conditions. These are included in the formulation using a penalty method, 
as in a standard Finite Element Method. 

Finally, the linear system of Eq. (41) is solved with a standard LU decomposition. The 
velocities computed in the {Utot} are used to update the grains positions at the next time 
step. 



3.6 Detection of contact 

The present model is intended to be used for the study of hot cracking, i.e., for Qg > 0.9, 
so that the width of the liquid channels is very small compared to the grain size, and the 
displacement of the grains will be very limited. It is important in the implementation of 
the model to detect and account for contact between grains. Equation (9) shows that the 
pressure drop in a channel tends to infinity as (l//i)^ as /i — 0, and this is sufficient to prevent 
the interpenetration of the grains in the simulation, as long as the time step is sufficiently 
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small. In some highly constrained situations, we found it necessary to implement a dynamic 
time step refinement procedure in the code in order to prevent grain interpenetration. 

It is more problematic to handle the changes of the environment of the grains, i.e., 
modifications of their nearest neighbors. This typically occurs after large displacements of 
the grains, especially for channels that are short to begin with. There are a few additional 
situations where the grains can interpenetrate (see [8] for further details). These cases are 
handled by revising the list of the first neighbors after each time step. [20]. In the example 
problems presented in the next section, the volume fraction of solid is high, and the system 
is thus sufficiently constrained that this phenomenon does not occur. 

4 Results 

4.1 Boundary conditions 

To use the model to investigate the mechanical behavior of the mushy zone, two test cases 
were analyzed where a tensile load was applied to a small ensemble of grains comprising 
a square volume element of edge Ly- Boundary conditions for the two cases are shown in 
Fig. 4. Since two phases are present, boundary conditions are needed for either the pressure 
or the fiow rate for the liquid channels, and for either the velocity or the applied forces on 
the solid grains. 



Ps=Po R=Po Ps =Po (I>r0 




(a)Ideal feeding (b)No feedin] 



Figure 4: Boundary conditions for two traction tests along the X axis: (a) Ideal feeding 
from the upper boundary; (b) No feeding. The grains or clusters of grains are colored with 
various grey levels. 

The two cases differ only in the boundary conditions applied on the top surface. In the 
first case (Fig. 4(a)), a pressure Pq is imposed on both the liquid and the solid, corresponding 
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to an average hydrostatic pressure on the system. ^ The fluid flow is free, thus allowing 
feeding from the upper boundary. We refer to this case as "ideal feeding" because it simulates 
a sample in contact with a liquid reservoir (a feeder) at pressure Pq. In the second case, we 
also impose a pressure Pq on the upper surface of the solid grains but the fluid flow is set to 
zero. This simulates a situation where feeding is impossible. 

The remaining boundary conditions are the same for both cases. The fluid flow and 
the horizontal X-component of the grain velocity are zero on the left boundary, while the 
y-component of the grain velocity is free, i.e., no forces along the vertical axis. This set of 
boundary conditions is equivalent to a symmetry plane. On the bottom boundary, the flux 
and the y-component of the velocity arc zero, while the X-component is free. On the right 
boundary, a velocity eLy is imposed on the solid in the X-direction to study the effect of an 
imposed strain rate e. The fluid flow is zero on the right boundary. ^ 

We note that these boundary conditions are similar to those of the model derived by 
Lahaie and Bouchard for a regular arrangement of hexagonal grains [35]. In our numerical 
calculations, the solidification of the system is calculated first, and then the mechanics of 
the mushy zone is computed without allowing any further solidification, i.e., at fixed solid 
fraction. This implies that gs is fixed, so that there is no solidification shrinkage. 

4.2 Tension tests 

Figure 5(a) shows the stress-strain curves for a sample solidified at constant T = — 1 K s^^, 
up to three volume fractions of solid {qs = 0.92, Qs = 0.94, and Qs = 0.96). Each sample 
was then strained along the X-direction at a rate e = 4 x 10~^ s~^. The results for the 
two tests, ideal and no feeding, are shown with open and filled symbols, respectively As 
can be seen, the two tests give the same stress-strain response at strains up to about 2.5% 
for Qg = 0.92, 1.2% for gg = 0.94 and 0.3% for gg = 0.96. Beyond these strains, the stress 
increases abruptly when feeding is not allowed, whereas it remains low and even decreases 
past a maximum in the case of ideal feeding. 

These results can be understood by looking at the local deformation mechanisms shown 
in Fig. 6, corresponding to the isothermal mushy zone strained at gg — 0.92. This domain 
contains 200 grains having an average diameter of 100 iim. Note that the computation of a 
traction test simulation on such a mushy zone takes about 30 s on a personal computer with 
a 2 GHz Intel Core Duo processor. Figure 6(1) shows the grain structure with the liquid 

^To represent the hydrostatic pressure on the solid, a force is imposed on each grain at the boundary. 
This force is oriented along the normal to the boundary and is equal to —PqLi, where Lh is the length of the 
external boundary of the grain. 

■^Note that for boundary conditions, the liquid flux is considered in the frame of the solid grains. In the 

laboratory frame, a fluid flux is observed due to the advection of the solid. 
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(a) (b) 

Figure 5: Stress as a function of strain for a square domain strained at a rate e — A x 10~^ 
s~^. Open and filled symbols correspond to ideal and no feeding, respectively, (a) Qs < 0.97, 
(b) Qs > 0.97 . 



channels at the onset of deformation, i.e., for e = 0. The grains are identified by assigning 
various gray levels, and the velocity of each grain is displayed with small arrows in Fig. 6(la) 
(the scale at the bottom of the figure gives the modulus of the velocity). When two grains 
establish a solid bond, i.e., the width of the corresponding channel goes to zero, they are 
shaded with the same gray scale, making it easy to recognize the formation of grain clusters 
as sohdification progresses. [25] At this relatively low volume fraction of sohd {qs — 0.92), 
only a few such clusters of grains have formed. 

In Fig. 6 (lb), the fluid flow in each channel is represented by lines whose width is pro- 
portional to their magnitude, the direction being indicated by a black triangle. The scale 
used for this flow representation is again shown at the bottom of the figure. As two grains 
get closer together, they squeeze the liquid out of the channel, and on the other hand, as 
they move away from each other, liquid is pumped into the channel. In both cases, the width 
of the corresponding line varies along the channel length. The flow can also be important 
in channels of flxed width that feed other regions of the mushy zone. As can be seen in 
Fig. 6 (lb), deformation of the mush is accommodated by fluid flow, but these flows remain 
small and localized. Long-range feeding of the mush is not necessary at low strain, and 
therefore both cases give the same behavior of the mush (see Fig. 5(a)). 

The stresses in the sample are also very low at low strain, as indicated in Fig. 6(lc). 
The interaction forces between the grains via the liquid channel are represented by a line 
connecting nucleation centers, whose width is proportional to the magnitude of the force 
(scale shown at the bottom), and whose gray scale level (or colour) indicates whether the 



16 



(3) 




(a)ciusters structure and velocity (b) Localisation of fluid flow (c) Forces between grains 



Figure 6: Isothermal mushy zone with = 0.92 and deformed at a strain rate e = 4 x 10^^ 
s~^ along the horizontal X-direction: (1) Mushy zone at 0% strain; (2) Mushy zone at 4% 
strain with ideal feeding from the upper face; (3) Mushy zone at 4% strain without feeding 
from the upper face. In (a), the grains (or grain clusters) are shown with various grey levels 
together with their velocity represented with small arrows. In (b), the flow in the channels 
is represented with a line of thickness proportional to the intensity and a triangle indicating 
the direction. In (c), the forces between the grains are represented with lines of variable 
thickness proportional to the modulus. A grey hue (red) corresponds to traction and a dark 
hue (blue) to compression. 
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force is tensile or compressive. At low volume fraction of solid and strain, the forces are of 
course very low and the corresponding lines are barely visible in Fig. 6(lc). 

Figures 6(2) shows the same mushy zone after 4% deformation for the ideal feeding case, 
and Fig. 6(3) shows the case where feeding is prohibited. The overall deformation is indicated 
in (a), the initial volume element is drawn as a dashed-line square. It is important to note 
that deformation is localized to a few channels, roughly oriented normally to the tensile 
direction. The fluid tends to flow from channels oriented in the direction of the stress to 
channels oriented perpendicular to the stress. The grains tend to be pulled inward along 
the vertical y-direction as we try to pull the mushy zone in the X-direction. This is also 
reflected by the forces shown in Fig. 6 (2c) and (3c), which essentially correspond to traction 
along the horizontal X-direction and to compression along the F-direction. 

It is also interesting to note that fluid flow is much more important in the case of ideal 
feeding (Fig. 6(2b)) compared to the case of no feeding (Fig. 6(3b)), even though the imposed 
strain rate is the same. This shows that redistribution of fluid occurs over larger distances 
with the accumulation of deformation. Fluid flow from the upper boundary is clearly visible 
in Fig. 6 (2b), and this flow relaxes the stresses in the upper part of the sample (Fig. 6 (2c)). 

If the mush cannot be fed from the upper boundary (Fig. 6(3)), redistribution of the fluid 
channels and of the grains also occurs as deformation increases. However, since no feeding 
from the top surface is allowed, the fluid follows a more difficult path, and the stresses in 
the samples are higher. Note that, as no hquid flow is allowed on any boundary for this case 
and the solid grains are rigid, the total volume of the specimen (grains + liquid) is constant. 

In summary, at low strain, deformation is accommodated by local redistribution of the 
liquid. This deformation is localized in liquid channels oriented roughly perpendicular to the 
stress direction. As deformation increases, more channels get closed and fluid redistribution 
occurs at a larger scale. At that point, the ability to feed the mush from some liquid reservoir 
becomes important. 

4.3 Transition in the feeding mechanism 

As the volume fraction of solid in the volume element increases, the difference between the 
ideal-feeding and the no-feeding cases occurs at lower strain (Fig. 5(a)). The reason for this 
is fairly obvious, as the width of the liquid channels decreases with increasing gg. 

However, for Qs > 0.97 (Fig. 5(b)), the behavior of the two cases becomes different even 
at the limit £ — > 0. To understand this new response, consider Fig. 7, which shows a 
mushy zone with Qg — 0.975 deformed at the same strain rate e — A x 10~^ s~^. At the 
onset of deformation (e — 0, Fig. 7(la)), it can be seen that the sohd grains start to form 
significant numbers of grain clusters. Moreover, Fig. 7(lb) shows that there is a significant 
redistribution of fiuid, even though the deformation of the mush is zero at that stage. The 
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fluid is located predominantly at the boundaries of the clusters which move as a single larger 
grain. It is interesting to compare this figure with Fig. 6(lb), which represents a mushy zone 
under the same straining and feeding conditions, but at a lower solid fraction. 




Clusters structure and velocity 



Localisation of fluid flow 



Figure 7: Isothermal mushy zone with Qs = 0.975 strained with e = 4 x 10~^ s~^. Feeding 
from the upper face is allowed: e = (1), e = 0.05% (2). In (a), the grains and grain clusters 
are shown with various grey levels and their velocity is indicated with small arrows. In (b), 
the flow in the channels is displayed with lines of thickness proportional to the intensity and 
triangles indicating the direction. 



The same mushy zone is represented in Fig. 7(2) at 0.05% deformation , i.e., after 0.125 
s. The very small deformation is localized to just a few channels that are located at the 
edges of the clusters and oriented normal to the tensile axis. Local redistribution of the 
liquid is difficult and most of the fluid to accommodate deformation is brought from the 
upper boundary even at this low strain level. Therefore, above Qs = 0.97, the deformation 
mechanism due to local fluid redistribution is no longer possible, due to the presence of 
relatively large clusters and thinner liquid channels. It is interesting to note that the volume 
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fraction of solid at which the mechanical behavior of the mush changes corresponds precisely 
to the value gs,i%iic at which 1% of the liquid channels become isolated and the permeability 
of the mush deviates from the Kozeny-Carman relationship. [25] 

The mechanical model we have presented shows that this solid fraction also corresponds 
to the point where accommodation of deformation by fluid flow becomes extremely difficult. 
Therefore, this point can be associated with the ductility minimum point observed experi- 
mentally [6,13], also called coalescence sohd fraction by some authors [7]. We find this point 
at a value Qs = 0.97, whereas the experiments performed on inoculated globular microstruc- 
turcs give instead a value Qs = 0.95 [7]. This difference is probably due to the 2D nature of 
our model, which tends to underestimate the size of the last liquid films (segments in a 2D 
Voronoi model instead of polyhedral surfaces in 3D). 

4.4 Discussion 

The effect of external loads on the response in the model can be deduced easily from the 
linearity of the matrix [E] (Eq. (29)). In particular, an increase of the metallostatic pressure 
Pq applied on the upper surface of Fig. 4 simply shifts the stress response of the solid 
network.^ This response corresponds to the Terzaghi effective stress, i.e., the behavior of 
the mush depends only on the difference between the apphed stress and the hydrostatic 
pressure [36] . Similarly, the stress response of the mushy zone varies linearly with the strain 
rate. 

On the other hand, the evolution of the mushy zone with strain requires a numerical cal- 
culation. The stress-strain curves in Fig. 5 can be compared qualitatively with experimental 
data (see e.g. [6-8]). For solid fractions lower than gs,i%iic (Fig. 5(a)), the strain at which 
the stress increases abruptly corresponds well to the experimental strains at fracture for the 
corresponding value of gg- However, the shape of the stress-strain curve and the magnitude 
of the stress are clearly different from the experimental values. At higher solid fractions, the 
strain at which the stress increases abruptly no longer correlates well with the experimental 
observations. 

We have not considered solid deformation in the model. For solid fraction lower than 
the ductility minimum {gs < gs,i%iic), solid deformation certainly plays a role in the overall 
behavior of the mush, but grain displacement is the dominant deformation mechanism. The 
model reproduces the strain at which interlocking of the grains occurs, even though rota- 
tions were not considered, but it cannot reproduce the shape of the stress-strain curve. At 
higher solid fractions {gg > gs,i%iic), solid deformation is clearly the dominant mechanism. 

^Prom a numerical point of view, it is clear that zero grain displacement and a uniform pressure is a 
solution of the problem. Because the system is linear, this implies that any solution can be shifted by a 
uniform amount of the pressure without affecting the displacement of the grains. 
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The ductility increase observed experimentally at high gg-valne corresponds to the strength 
increase of the solid network due to its progressive percolation. 

This 2D model demonstrates a transition in the mechanical behaviour of the mush as the 
solid fraction increases. We have included the 3D aspects of that transition in a simplified 
way. The details of this transition requires a fully 3D model, which we will describe in a 
forthcoming publication. 

4.5 First consequences on hot tearing criteria 

Hot tearing criteria based on the feeding ability of the mush [11, 12, 16] predict the hot 
cracking sensitivity (HCS) of alloys fairly well. Such criteria are based on a critical strain rate, 
which seems to be particularly applicable in processes such as DC casting, where thermally- 
induced stresses are essentially perpendicular to the thermal gradient. [37,38] However, in 
tensile test experiments where the applied stress is parallel to the thermal gradient, hot 
cracking is found to be largely independent of the stain rate. [8, 38] As mentioned in the 
previous section, the response of the present model is at first sight proportional to the strain 
rate. However, the granular nature of the model makes it strongly strain-dependent. 



Gas 




Gas - t -] Liquid 


(a) Solid 


(b) Solid 



Figure 8: (a) Equilibrium of forces at triple junction, (b) Shape of the meniscus between 
two grains. 



The criterion for hot tear nucleation is usually written in the form 

Pi < Pc (42) 

where pi is the local pressure in the liquid and Pc represents the cavitation pressure at which 
a pore nucleates. As a first approximation, this cavitation pressure can be estimated as the 
overpressure required to overcome capillary forces at the liquid-pore interface. Since pores 
will form at high solid fraction, i.e., in very narrow liquid channels, their radius of curvature 
is dictated by the width of the channel, as illustrated in Fig. 8. At the triple junction between 
the liquid, solid and pore, the equilibrium condition is given by: 

7g/ COS e = 'ygs- J si (43) 
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Figure 9: Maximum value of the cavitation pressure in the mushy zone as a function of strain 
for an impregnation factor / of 1 Jm~^. 



where 7^^, 7^; and 7^^ are the interfacial energy between sohd and hquid, hquid and pore, 
and sohd and pore, respectively, and G is the dihedral angle. Therefore, the radius of the 
pore R is given by: 

(44) 

cos © 

where h is the half width of the hquid channel (see Fig. 8 (b)). Therefore 

where / is called the "impregnation factor." Since liquid metals wet very well their own 
solid, i.e., 'jsi <^ Igi, the value of / is typically close to 7^/, e.g., 1 J m~^ for Al. 

Using such relationships, the maximal cavitation pressure Pc,max is plotted in Fig. 9 as 
a function of deformation for isothermal mushy zones with various solid fractions deformed 
at a strain rate £ = 4 x 10^'^ s^^. The pressure corresponds to the pressure necessary for 
cavitation of a pore in the widest channel of the mush, hmax- This channel width decreases 
with increasing Qs and considerably increases with strain, in particular for high solid fraction 
samples. Note that this evolution is largely independent of the other parameters. 

The variations of Pc due to strain are therefore more important than those induced by 
the strain rate. This phenomenon can explain the apparent insensitivity of hot cracking to 
strain rate, at least for tensile test experiments. This idea is further illustrated in Fig. 10 
where a mushy zone is represented at various strain levels. In this test, the strain rate i. 
is equal to 4 x 10"^ s~^, Qs — 0.92 and feeding from the upper boundary is not allowed. 
The localization of the fluid flow and the grain velocity are represented on the same picture. 
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Channels in which feeding is not represented (white channels) , correspond to those in which a 
pore has nucleated. In order to reach depressions capable of producing a pore by cavitation, 
the pressure at the upper boundary was fixed to Pq = —120 kPa. The impregnation factor 
I in these simulations was fixed to 1 J m~^. 

Figure 10 clearly shows the nucleation of pores in the largest channels, where deformation 
has been localized. Note that this approach allows an estimation of the appearance of damage 
in the mushy zone, but cannot model fracture, which would require explicit modeling of the 
deformation of the solid grains. 




Figure 10: Mushy zone at various strain levels with Qs = 0.92 and e = A x 10^^ s~^. Feeding 
from the upper face is not allowed. Grain velocity (black arrows) and fluid flow (grey/red 
lines) are represented on the same picture. The channels on which feeding is not represented 
(white channels) correspond to those where a pore has nucleated. 



5 Conclusion 

This paper has presented the development of a 2D granular model for the mechanical behavior 
of mushy zones that brings valuable insight into the interrelations between intergranular flow 
and grain movement. The model is based on conservation of mass and force at the level of 
the intergranular region. It produces naturally a form that is consistent with non-equilibrium 
thermodynamics. The whole mechanical problem can thus be expressed as a minimization 
of dissipation. Therefore, the results of the present simulations are well suited for continuum 
models based on a dissipation potential. [7] 

The model accounts directly for the random nature of grain nucleation, and for the 
progressive formation of grain clusters during solidification. [25] In addition, the granular 
nature of the material is modeled, which leads to further localization of deformation and 
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feeding (besides the localization due to grain clusters) . Both of these important phenomena 
cannot be resolved using typical volume-averaged formulations that smear out the details of 
the grain structure. 

The formulation of the model requires only a few integration points for each grain. Com- 
putation times are therefore very low and leave room for further extensions, in particular to 
three dimensions. 
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A Scaling analysis and simplification of the equations 

The X— component of the momentum balance for a constant density Newtonian fluid in a 
liquid channel is given by 

[dvx , dvx , dvx\ Dvx dp d^vx , d'^vx. 

We follow the procedures in [39] to put this equation in dimensionless form. To that end, 
we define the following scaled variables 

^ ^z;' ^ =2^' = — ^^ = -17^' ' -t:^p = ap7 ^^^^ 

where = {V^ - V^), V = {V^ - V^) and AP^ = {Pj - Pi)- Substituting these scaled 
variables into Eq. (46) gives 

IxL^ Dt° ixVtL^ dX° LI dX°'^ dY°'^ ^ ^ 

The factor 2h/Lc that appears in several terms can be expressed in terms of the solid fraction 

^tot \ gs' J 

where Htot is the height of the elementary triangle and Ltot the length of its base. For a 
regular hexagonal network of solid grains this is constant, given by 

which is precisely the term introduced by Lahaie and Bouchard in their hot tearing criterion 

_ 1 

[35]. We are interested in solid fractions Qs > 0.8 for which {qs ^ — 1) ~ (1 — 5's)/2 < 0.1. 
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To estimate the relative magnitude of the various terms, we consider the typical values 
of the physical parameters in an inoculated Al-Cu alloy pi ~ 2440 Kg m~^, // ~ 1.5 x 10~^ 
Pa s, Lc ~ 10~^ m. [40] Moreover, in DC casting the typical strain rate is on the order of 
e ~ 10-3 s-^ [2] and thus V* ~ L^e ~ 10"^ m s^^ Thus 

^ ~ 10- « 1 (51) 

and the transient and inertial terms of Eq. (48) are negligible. Similarly: 

— ~ 10-2 < 1 52 

and we can reasonably neglect the term in d'^Vx/dX'^. This leaves the simplified equation 
describing the flow in a channel: 

Applying the same procedure to the equation for the y-component of the fluid velocity gives: 

d'^v^ _ 2/iAPy dp° 

To estimate of the pressure variation in the channel, we set the dimensionless groups to 
one [39] 

and thus 

Thus, to a flrst approximation, we can neglect the variation of pressure across the channel in 
comparison to the variation of pressure along the channel. Therefore, even though fluid flow 
exists in the y-direction due to solidiflcation shrinkage or grain displacement, the Poiseuille 
equation gives a good approximation of the flow. 

In the derivation of the model, we considered only the dissipation along the channels 
and neglected the dissipation at triple junctions. Indeed, the dissipation along the channels 
is mainly due to viscous losses, whereas the dissipation at triple junctions is due to the 
direction change of the flow. Such losses are proportional to the kinetic energy of the flow 
[41]. However, Eq. (51) shows that the ratio between kinetic energy and viscous dissipation 
(Reynolds number) is on the order of 10-^ and therefore it is reasonable to neglect the 
dissipation at triple junctions. 
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B Detailed integration of the constitutive equations 



B.l Basic equations 

Considering that dp/dX is constant along the width of the channel (see Appendix A), the 
integration of Eq. (3) (see main section) gives: 

where a non-slip condition is considered at the s — I interfaces. Therefore, we have: 



and from Eq. (6) we get: 



2h^ 



^ (59) 



3/i dX^ 

where the source term V"' for the flow in the channel is given by: 

= - V: = VS, - VS^ + 2Pv* (60) 
Finally, the pressure and velocity profiles are given by 

P{X) = ^(X^ - (^)^) + + (61) 



The integration of Eq. (62) gives the fluid flux in the channel: 

$,_,(X) = -V-X +^{P^- Pj) + 2h^^ (63) 

B.2 Liquid flux 

The mass balance in each channel must be completed by a mass balance at each vertex. As 
the fluid is considered as incompressible, one has 

s 

where the summation is carried out over the three channels s = 1, 2, 3 connected to vertex i 
(see Fig. 3). It is possible to show that this summation is equivalent to [8] 

E |^(^^ - + 2/^^^^^^ + - L^aV: = (65) 

where the indices of the grains (a and h) correspond to the neighbors of each channel. 

26 



B.3 Forces 



In order to get a representation of the stress tensor, let us consider the pressure variation in 
the F-direction. For an incompressible fluid, we have 

V • {v} = (66) 

As the speed profile along the X-axis is already defined (Eq. (62)), we get 

dvY _ _dvx _ _SV" 
~dY ~ ~~dX ~ ~4/^ 

Thus 



{Y' - h') (67) 



The variation of pressure in the y-direction is obtained from the momentum balance equa- 
tion, 

dp d'^VY 

dY^^W^ (69) 
We get the variation of pressure along the F-direction 

oT/n 

P = + (70) 

where f{X) is the pressure in the channel for the line y = which we take equal to the 
pressure profile of Eq. (61). Finally the pressure profile is: 

p(x.r)^^fA-^-f^r-yA.^x.^ (71, 



4/i3 I \2 

As shown in the scaling analysis, the pressure loss in the F-direction is negligible in front 
of the pressure variation in the X-direction. Note that we found here the relationship 
APy/APx ~ /Li whereas the scahng analysis gives a relationship APy/APx ~ h/Lc 
(Eq. 56). 

The stress tensor in the liquid can be derived from 



-p+2^l^ Mfr + ff) 



('^' ^) I / dvx _i_ dvY \ „ _i_ o , , dvY ' (''^^) 



and thus 



-P-'-i^(Y^-l' 



where p is given by Eq. (71). Note that the dxvy term is nil. This stress field is coherent as 
we can check that : 

V • [a] (X, Y) = {0} (74) 
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For the simplicity of the equations, we choose to neglect the variation of pressure in the 
y-direction (Eq. (61)). Thus, we should also neglect the Y"^ terms in the diagonal of the 
stress tensor as they are exactly of the same order. 



-P 

Note that the divergence of this tensor is not zero. However, considering this stress tensor, 
the sum of forces and of rotational momentum on the boundaries of a channel are zero, which 
is the condition for the coherency of the numerical scheme. 

Equation (75) gives the stress tensor in the part of the channel where the two grains 
match. In the part where they do not match, we simply consider that we have a homogeneous 
pressure equal to the pressure of the integration point (Pj or Pj), see Fig. 2). Thus, we can 
integrate the force exerted by a grain on the liquid 



and 



,F , _ I (P'- P'V' + "ItW - : ,77, 
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